  subroutine radar_simulator( &
    hp, &								! hydrometeor and radar simulator properties 
    nprof,ngate, &							! size of input matricies
    hgt_matrix,hm_matrix,re_matrix,Np_matrix, &				! inputs
    p_matrix,t_matrix,rh_matrix, &					
    Ze_non,Ze_ray,h_atten_to_vol,g_atten_to_vol,dBZe, &			! outputs
    g_to_vol_in,g_to_vol_out)						! optional in/out
 
  use m_mrgrnk 
  use array_lib
  use math_lib
  use optics_lib
  use radar_simulator_types
  use scale_LUTs_io
  implicit none
  
! Purpose:
!
!   Simulates a vertical profile of radar reflectivity
!   Originally Part of QuickBeam v1.04 by John Haynes & Roger Marchand.
!   but has been substantially modified since that time by
!   Laura Fowler and Roger Marchand (see modifications below).
!
! Inputs:
!
!   [hp]              structure that defines hydrometeor types and other radar properties
!
!   [nprof]           number of hydrometeor profiles
!   [ngate]           number of vertical layers
!
!   (The following 5 arrays must be in order from closest to the radar
!    to farthest...)
!
!   [hgt_matrix]      height of hydrometeors (km)
!   [p_matrix]        pressure profile (hPa)
!   [t_matrix]        temperature profile (C)
!   [rh_matrix]       relative humidity profile (%) -- only needed if gaseous aborption calculated.
!
!   [hm_matrix]       table of hydrometeor mixing rations (g/kg)
!   [re_matrix]       OPTIONAL table of hydrometeor effective radii.  0 ==> use defaults. (units=microns)    
!   [Np_matrix]	      Optional table of hydrometeor number concentration.  0 ==> use defaults. (units = 1/kg)
!
! Outputs:
!
!   [Ze_non]          radar reflectivity without attenuation (dBZ)
!   [Ze_ray]          Rayleigh reflectivity (dBZ)
!   [h_atten_to_vol]  attenuation by hydromets, radar to vol (dB)
!   [g_atten_to_vol]  gaseous atteunation, radar to vol (dB)
!   [dBZe]            effective radar reflectivity factor (dBZ)
!
! Optional:
!   [g_to_vol_in]     integrated atten due to gases, r>v (dB).
!                     If present then is used as gaseous absorption, independently of the
!                     value in use_gas_abs
!   [g_to_vol_out]    integrated atten due to gases, r>v (dB).
!                     If present then gaseous absorption for each profile is returned here.
!
! Created:
!   11/28/2005  John Haynes (haynes@atmos.colostate.edu)
!
! Modified:  
!   09/2006  placed into subroutine form (Roger Marchand,JMH)
!   08/2007  added equivalent volume spheres, Z and N scalling most distrubtion types (Roger Marchand)
!   01/2008  'Do while' to determine if hydrometeor(s) present in volume
!             changed for vectorization purposes (A. Bodas-Salcedo)
!
!   07/2010  V3.0 ... Modified to load or save scale factors to disk as a Look-Up Table (LUT)
!		  ... All hydrometeor and radar simulator properties now included in hp structure
!		  ... hp structure should be initialized by call to radar_simulator_init prior 
!		  ... to calling this subroutine.  
! 	     Also ... Support of Morrison 2-moment style microphyscis (Np_matrix) added 
!  		  ... Changes implement by Roj Marchand following work by Laura Fowler
!
!   10/2011  Modified ngate loop to go in either direction depending on flag 
!	     hp%radar_at_layer_one.  This affects the direction in which attenuation is summed.
!
!	     Also removed called to AVINT for gas and hydrometeor attenuation and replaced with simple
!	     summation. (Roger Marchand)
!
!		
! ----- INPUTS -----  
 
  logical, parameter  ::  DO_LUT_TEST = .false. 
  logical, parameter  ::  DO_NP_TEST = .false. 

  type(class_param), intent(inout) :: hp
  
  integer, intent(in) ::  nprof,ngate
    
  real*8, dimension(nprof,ngate), intent(in) :: &
  	hgt_matrix, p_matrix,t_matrix,rh_matrix
  	
  real*8, dimension(hp%nhclass,nprof,ngate), intent(in) :: hm_matrix
  real*8, dimension(hp%nhclass,nprof,ngate), intent(inout) :: re_matrix
  real*8, optional, dimension(hp%nhclass,nprof,ngate), intent(in)    :: Np_matrix
    
! ----- OUTPUTS -----
  real*8, dimension(nprof,ngate), intent(out) :: Ze_non,Ze_ray, &
 	g_atten_to_vol,dBZe,h_atten_to_vol

! ----- OPTIONAL -----
  real*8, optional, dimension(ngate,nprof) :: &
  g_to_vol_in,g_to_vol_out ! integrated atten due to gases, r>v (dB). This allows to output and then input
                           ! the same gaseous absorption in different calls. Optional to allow compatibility
                           ! with original version. A. Bodas April 2008.
        
!  real*8, dimension(nprof,ngate) :: kr_matrix 

! ----- INTERNAL -----

  integer :: &
  phase, &			! 0=liquid, 1=ice
  ns 				! number of discrete drop sizes

  integer*4, dimension(ngate) :: &
  hydro				! 1=hydrometeor in vol, 0=none
  real*8 :: &
  rho_a, &			! air density (kg m^-3)
  gases				! function: 2-way gas atten (dB/km)

  real*8, dimension(:), allocatable :: &
  Di, Deq, &      		! discrete drop sizes (um)
  Ni, &    			! discrete concentrations (cm^-3 um^-1)
  rhoi				! discrete densities (kg m^-3)
  
  real*8, dimension(ngate) :: &
  z_vol, &			! effective reflectivity factor (mm^6/m^3)
  z_ray, &                      ! reflectivity factor, Rayleigh only (mm^6/m^3)
  kr_vol, &			! attenuation coefficient hydro (dB/km)
  g_vol, &			! attenuation coefficient gases (dB/km)
  a_to_vol, &			! integrated atten due to hydometeors, r>v (dB)
  g_to_vol			! integrated atten due to gases, r>v (dB)
   
 
  integer,parameter :: KR8 = selected_real_kind(15,300)
  real*8, parameter :: xx = -1.0_KR8
  real*8,  dimension(:), allocatable :: xxa
  real*8 :: kr, ze, zr, pi, scale_factor, tc, Re, ld, tmp1, ze2, kr2, apm, bpm
  real*8 :: half_a_atten_current,half_a_atten_above
  real*8 :: half_g_atten_current,half_g_atten_above
  integer*4 :: tp, i, j, k, pr, itt, iff

  real*8    step,base, Np
  integer*4 iRe_type,n,max_bin
 
  integer   start_gate,end_gate,d_gate
  
  logical :: g_to_vol_in_present, g_to_vol_out_present
	
  ! Logicals to avoid calling present within the loops
  g_to_vol_in_present  = present(g_to_vol_in)
  g_to_vol_out_present = present(g_to_vol_out)
  
  !
  ! load scaling matricies from disk -- but only the first time this subroutine is called
  !
  if(hp%load_scale_LUTs) then
  	 	
  	call load_scale_LUTs(hp)
  	hp%load_scale_LUTs=.false.  	
   	hp%Z_scale_added_flag = .false.		! will be set true if scaling Look Up Tables are modified during run	
  
  endif
  
  pi = acos(-1.0)
	
  ! // loop over each profile (nprof)
  do pr=1,nprof

!   ----- calculations for each volume ----- 
    z_vol(:) = 0
    z_ray(:) = 0
    kr_vol(:) = 0
    hydro(:) = 0    

!   // loop over each range gate (ngate) ... starting with layer closest to the radar !

    if(hp%radar_at_layer_one) then
    	start_gate=1
    	end_gate=ngate
    	d_gate=1
    else
      	start_gate=ngate
    	end_gate=1
    	d_gate=-1
    endif

    do k=start_gate,end_gate,d_gate
    
!     :: determine if hydrometeor(s) present in volume
      hydro(k) = 0
      do j=1,hp%nhclass ! Do while changed for vectorization purposes (A. B-S)
        if ((hm_matrix(j,pr,k) > 1E-12) .and. (hp%dtype(j) > 0)) then
          hydro(k) = 1
          exit
        endif
      enddo

!     :: if there is hydrometeor in the volume            
      if (hydro(k) == 1) then

        rho_a = (p_matrix(pr,k)*100.)/(287*(t_matrix(pr,k)+273.15))

	if(rho_a<0 .or. rho_a>2) then
	
		write(*,*) pr,k,p_matrix(pr,k),t_matrix(pr,k)
		write(*,*) 
	endif

!       :: loop over hydrometeor type
        do tp=1,hp%nhclass

          if (hm_matrix(tp,pr,k) <= 1E-12) cycle

	  phase = hp%phase(tp)
	  if(phase==0) then
		itt = infind(hp%mt_ttl,t_matrix(pr,k))
  	  else
		itt = infind(hp%mt_tti,t_matrix(pr,k))
          endif

  	  if(re_matrix(tp,pr,k).eq.0) then
			
	  	call calc_Re(hm_matrix(tp,pr,k),Np_matrix(tp,pr,k),rho_a, &
  		     	hp%dtype(tp), &
  		     	hp%dmin(tp),hp%dmax(tp), &
  		     	hp%apm(tp),hp%bpm(tp), &
  		     	hp%rho(tp),hp%p1(tp),hp%p2(tp),hp%p3(tp), &
  		     Re)
  		     
  		     ! if(tp.eq.9) then
  		     !		write(*,*) tp,pr,k,hm_matrix(tp,pr,k),Np_matrix(tp,pr,k),Re
  		     ! endif
 
		re_matrix(tp,pr,k)=Re
	  else
	  	if(Np_matrix(tp,pr,k)>0) then
	  		print *, 'Warning: Re and Np set for the same ', &
	  			 'volume & hydrometeor type.  Np is being ignored.'
	  	endif
	  	
	  	Re = re_matrix(tp,pr,k)
	  endif	
	    
	  iRe_type=1
	  
	  if(Re.gt.0) then

	  	! determine index in to scale LUT
		!
		! distance between Re points (defined by "base" and "step") for
		! each interval of size Re_BIN_LENGTH
		
		n=floor(Re/Re_BIN_LENGTH)
		
		if(n>=Re_MAX_BIN) then
			n=Re_MAX_BIN-1	
		endif
		step=hp%step_list(n+1)
		base=hp%base_list(n+1)
		
		iRe_type=floor(Re/step)

		if(iRe_type.lt.1) then  
			iRe_type=1			
		endif

		Re=step*(iRe_type+0.5)	! set value of Re to closest value in allowed in LUT.
		iRe_type=iRe_type+base-floor(n*Re_BIN_LENGTH/step)
	
	 	! make sure iRe_type is within bounds
		if(iRe_type.ge.nRe_types) then  

			write(*,*) 'Warning: size of Re exceed value permitted ', &
				   'in Look-Up Table (LUT).  Will calculate. ' ! , &
				   ! Re,iRe_type,nRe_types,step,base

			! no scaling allowed
			iRe_type=nRe_types
			hp%Z_scale_flag(tp,itt,iRe_type)=.false.
		else
			! set value in re_matrix to closest values in LUT
			if(.not. DO_LUT_TEST) then
				re_matrix(tp,pr,k)=Re
			endif
		
		endif
	  endif
	
  	  ! use Ze_scaled, Zr_scaled, and kr_scaled ... if know them
	  ! if not we will calculate Ze, Zr, and Kr from the distribution parameters
  	  if( (.not. hp%Z_scale_flag(tp,itt,iRe_type)) .or. DO_LUT_TEST)  then
 	 
	  	! :: create a distribution of hydrometeors within volume	  
	  	select case(hp%dtype(tp))
          	case(4)
		    ns = 1
		    allocate(Di(ns),Ni(ns),rhoi(ns),xxa(ns),Deq(ns))
		    Di = hp%p1(tp)
		    Ni = 0.
		case default
 		    ns = nd   ! constant defined in radar_simulator_types.f90            
		    allocate(Di(ns),Ni(ns),rhoi(ns),xxa(ns),Deq(ns))
 		    Di = hp%D
 		    Ni = 0.
		end select

	  	call dsd(hm_matrix(tp,pr,k),re_matrix(tp,pr,k),Np_matrix(tp,pr,k), &
	  		Di,Ni,ns, &
	  		hp%dtype(tp),rho_a,t_matrix(pr,k), &
	  		hp%dmin(tp),hp%dmax(tp),hp%apm(tp),hp%bpm(tp), &
	    		hp%rho(tp),hp%p1(tp),hp%p2(tp),hp%p3(tp))

		! calculate particle density 
	  	if (phase == 1) then
	  	
	    		if (hp%rho(tp) < 0) then
                
				! MG Mie approach - adjust density of sphere with D = D_characteristic to match particle density		
				! hp%rho_eff(tp,1:ns,iRe_type) = (6/pi)*hp%apm(tp)*(Di*1E-6)**(hp%bpm(tp)-3)   !MG Mie approach
		
				! as the particle size gets small it is possible that the mass to size relationship of 
				! (given by power law in hclass.data) can produce impossible results 
				! where the mass is larger than a solid sphere of ice.  
				! This loop ensures that no ice particle can have more mass/density larger than an ice sphere.
				! do i=1,ns
				! if(hp%rho_eff(tp,i,iRe_type) > 917 ) then
				!	hp%rho_eff(tp,i,iRe_type) = 917
				!endif
				!enddo

				! alternative is to use equivalent volume spheres.
			    	hp%rho_eff(tp,1:ns,iRe_type) = 917  				! solid ice == equivalent volume approach
			      	Deq = ( ( 6/pi*hp%apm(tp)/917 ) ** (1.0/3.0) ) * &
				   ( (Di*1E-6) ** (hp%bpm(tp)/3.0) )  * 1E6 		

            		else

        		    	! hp%rho_eff(tp,1:ns,iRe_type) = hp%rho(tp)   !MG Mie approach
	     	
				! Equivalent volume sphere (solid ice rho_ice=917 kg/m^3).
			     	hp%rho_eff(tp,1:ns,iRe_type) = 917
			     	Deq=Di * ((hp%rho(tp)/917)**(1.0/3.0))  

	    		endif
	    		
	    	else
	    		! I assume here that water phase dopplets are spheres.
	    		! hp%rho should be ~ 1000  or hp%apm=524 .and. hp%bpm=3
	    		Deq = Di;		
	  	endif
	  	  
		! calculate effective reflectivity factor of volume
	  	xxa = -9.9
	  	rhoi = hp%rho_eff(tp,1:ns,iRe_type)
	  	call zeff(hp%freq,Deq,Ni,ns,hp%k2,t_matrix(pr,k),phase,hp%do_ray, &
	      		ze,zr,kr,xxa,xxa,rhoi)
  	
  		! test code ... compare Np value input to routine with sum of DSD
  		! NOTE: if .not. DO_LUT_TEST, then you are checking the LUT approximation 
  		!       not just the DSD representation given by Ni
		if(Np_matrix(tp,pr,k)>0 .and. DO_NP_TEST ) then
			
			Np = path_integral(Ni,Di,1,ns-1)/rho_a*1E6
			
			! Note: Representation is not great or small Re < 2 
			if( (Np_matrix(tp,pr,k)-Np)/Np_matrix(tp,pr,k)>0.1 ) then !&
			    ! .and. 10*log10(ze)>-50 ) then
				
				write (*,*) 'Error: Np input does not match sum(N)'
			
				write(*,*) tp,pr,k,Re,Ni(1),Ni(ns),10*log10(ze)
				write(*,*) Np_matrix(tp,pr,k),Np, &
					   (Np_matrix(tp,pr,k)-Np)/Np_matrix(tp,pr,k)
				write(*,*)
			endif		
		endif
			
		deallocate(Di,Ni,rhoi,xxa,Deq)
  	
 	 	! LUT test code ... this segment of code compares full calculation to scaling result
  		if( hp%Z_scale_flag(tp,itt,iRe_type) .and. DO_LUT_TEST )  then
  	
	  		scale_factor=rho_a*hm_matrix(tp,pr,k) 

			! if more than 2 dBZe difference print error message/parameters.
			if( abs(10*log10(ze) - & 
			    10*log10(hp%Ze_scaled(tp,itt,iRe_type) * scale_factor)) &
			    > 2 ) then
			
				write(*,*) 'Roj Error: ',tp,itt,iRe_type,hp%Z_scale_flag(tp,itt,iRe_type),n,step,base
				write(*,*) 10*log10(ze),10*log10(hp%Ze_scaled(tp,itt,iRe_type) * scale_factor)
				write(*,*) hp%Ze_scaled(tp,itt,iRe_type),scale_factor
				write(*,*) re_matrix(tp,pr,k),Re
				write(*,*)				 
			endif
				
		endif
  	
	  else ! can use z scaling
	   
	    	scale_factor=rho_a*hm_matrix(tp,pr,k) 

		zr = hp%Zr_scaled(tp,itt,iRe_type) * scale_factor 
		ze = hp%Ze_scaled(tp,itt,iRe_type) * scale_factor
		kr = hp%kr_scaled(tp,itt,iRe_type) * scale_factor	

	  endif  ! end z_scaling
 
	  kr_vol(k) = kr_vol(k) + kr
	  z_vol(k) = z_vol(k) + ze
	  z_ray(k) = z_ray(k) + zr
	
	  ! construct Ze_scaled, Zr_scaled, and kr_scaled ... if we can
	  if( .not. hp%Z_scale_flag(tp,itt,iRe_type) ) then

		if(iRe_type>1) then
		
			scale_factor=rho_a*hm_matrix(tp,pr,k) 

			hp%Ze_scaled(tp,itt,iRe_type) = ze/ scale_factor
			hp%Zr_scaled(tp,itt,iRe_type) = zr/ scale_factor
			hp%kr_scaled(tp,itt,iRe_type) = kr/ scale_factor

			hp%Z_scale_flag(tp,itt,iRe_type)=.True.
			hp%Z_scale_added_flag(tp,itt,iRe_type)=.True.
		endif

	  endif

        enddo	! end loop of tp (hydrometeor type)

      else

	!     :: volume is hydrometeor-free	
 	kr_vol(k) = 0
	z_vol(k) = -999
        z_ray(k) = -999
	
      endif
      
       !     :: attenuation due to hydrometeors between radar and volume
       !
       if(d_gate==1) then
     	
		! a_to_vol(k) = 2*path_integral(kr_vol,hgt_matrix(pr,:),1,k-1)  ! Original integration ... Note requires level 1 be TOA (d_gate==1)
									        ! NOTE this integrages attenuation only for the layers ABOVE
									        ! the current layer ... i.e. 1 to k-1 rather than 1 to k ...
									        ! which may be a problem.   ROJ
									        ! in the new scheme I assign have the attenuation to the current layer
		
		! dheight cacluations assumes hgt_matrix points are the cell mid-points.
			 
     		if(k>2) then
	 		a_to_vol(k)=  a_to_vol(k-1) + &		! add to previous value to half of above layer + half of current layer
	 			      (kr_vol(k-1)+kr_vol(k))*(hgt_matrix(pr,k-1)-hgt_matrix(pr,k)) ; 	
	 		
	       	else
	  
	      	  	a_to_vol(k)=  kr_vol(k)*(hgt_matrix(pr,k)-hgt_matrix(pr,k+1));
	     	endif
   	   
   	     	
     	else   ! d_gate==-1
   
     		if(k<ngate) then
	 	     	a_to_vol(k) = a_to_vol(k+1) + & 	! add to previous value to half of above layer + half of current layer
	 			      (kr_vol(k+1)+kr_vol(k))*(hgt_matrix(pr,k+1)-hgt_matrix(pr,k)) ; 	
	       	    	
	       	else
	      		a_to_vol(k)= kr_vol(k)*(hgt_matrix(pr,k)-hgt_matrix(pr,k-1)); 
	     	endif	
	endif     
	
	 
      !     :: attenuation due to gaseous absorption between radar and volume
      if (g_to_vol_in_present) then
        g_to_vol(k) = g_to_vol_in(k,pr)
      else
        if ( (hp%use_gas_abs == 1) .or. ((hp%use_gas_abs == 2) .and. (pr == 1)) )  then
         
            g_vol(k) = gases(p_matrix(pr,k),t_matrix(pr,k)+273.15, &
            		     rh_matrix(pr,k),hp%freq)
            
            if(d_gate==1) then
     	
		! g_to_vol(k) = path_integral(g_vol,hgt_matrix(pr,:),1,k-1)  ! Original integration ... Note requires level 1 be TOA (d_gate==1)
     				 
     		if(k>1) then
	 		g_to_vol(k)=  g_to_vol(k-1) + &		! add to previous value to half of above layer + half of current layer
	 			      0.5*(g_vol(k-1)+g_vol(k))*(hgt_matrix(pr,k-1)-hgt_matrix(pr,k)) ; 	
	       	else
	      	  	g_to_vol(k)=  0.5*g_vol(k)*(hgt_matrix(pr,k)-hgt_matrix(pr,k+1));
	     	endif
     	    	
   	
     	    else   ! d_gate==-1
   
     		if(k<ngate) then
	 	     	g_to_vol(k) = g_to_vol(k+1) + & 	! add to previous value to half of above layer + half of current layer
	 			      0.5*(g_vol(k+1)+g_vol(k))*(hgt_matrix(pr,k+1)-hgt_matrix(pr,k)) ; 	
	       	    	
	       	else
	      		g_to_vol(k)= 0.5*g_vol(k)*(hgt_matrix(pr,k)-hgt_matrix(pr,k-1)); 
	     	endif	

	    endif     		
       
        elseif (hp%use_gas_abs == 0) then
            g_to_vol(k) = 0
        endif  
      endif
    	
      ! test code
      ! if(pr==1) then
      !  write(*,*) 'my test =',pr,k,kr_vol(k),a_to_vol(k),g_to_vol(k),hgt_matrix(pr,k) ! , &
       		! p_matrix(pr,k),t_matrix(pr,k)+273.15, &
                ! rh_matrix(pr,k),hp%freq
      !  endif            
	
      !     :: store results in matrix for return to calling program
      h_atten_to_vol(pr,k)=a_to_vol(k)
      g_atten_to_vol(pr,k)=g_to_vol(k)
      if ((hp%do_ray == 1) .and. (z_ray(k) > 0)) then
        Ze_ray(pr,k) = 10*log10(z_ray(k))
      else
        Ze_ray(pr,k) = -999
      endif
      if (z_vol(k) > 0) then
        dBZe(pr,k) = 10*log10(z_vol(k))-a_to_vol(k)-g_to_vol(k)
        Ze_non(pr,k) = 10*log10(z_vol(k))
      else
        dBZe(pr,k) = -999
        Ze_non(pr,k) = -999
      endif
      
    enddo	! end loop of k (range gate)
    
    ! Output array with gaseous absorption
    if (g_to_vol_out_present) g_to_vol_out(:,pr) = g_to_vol
    
  enddo		! end loop over pr (profile)  

  ! save any updates made 
  if(hp%update_scale_LUTs) then
	call save_scale_LUTs(hp)
  endif

  end subroutine radar_simulator
  
